About the Functional Form of the Parisi Overlap Distribution for the 
Three-Dimensional Edwards-Anderson Ising Spin Glass 
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Recently, it has been conjectured that the statistics of extremes is of relevance for a large class 
of correlated system. For certain probability densities this predicts the characteristic large x fall- 
off behavior f{x) ~ exp(— ae^), a > 0. Using a multicanonical Monte Carlo technique, we have 
calculated the Parisi overlap distribution P{q) for the three-dimensional Edward-Anderson Ising 
spin glass at and below the critical temperature, even where P{q) is exponentially small. We find 
that a probability distribution related to extreme order statistics gives an excellent description of 
P{q) over about 80 orders of magnitude. 

PACS: 75.10Nr,75.40Mg,75.50Lk. 
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The three-dimensional (3D) Edwards-Anderson |^] 
Ising (EAI) spin-glass model is a prototype of a disor- 
dered system, for which conflicting constraints create a 
rough free energy landscape. Such systems are of impor- 
tance for the understanding of a wide range of phenom- 
ena in physics, chemistry, biology and computer science. 
The overlap q between two replicas of the EAI model 
serves as an order parameter. Its probability density 
P{q) is, therefore, a quantity of central physical interest. 
More than twenty years ago, Parisi succeeded to calcu- 
late P{q) in the mean- field approximation However, 
for 3D physical systems the precise form of P{q) in the 
spin-glass phase, and the very nature of this phase, have 
remained a subject of debate More recently con- 

nections between spin glasses and extreme order statis- 
tics were pointed out by Bouchaud and Mezard . 
Then, Bramwell et al. ||] conjectured that the statistics 
of extremes, in particular Gumbel's first asymptote (in- 
troduced below), is of relevance for a large class of cor- 
related systems. We have used a multicanonical Monte 
Carlo (MC) technique to calculate P{q) numeri- 

cally, even when it is exponentially small. At the critical 
point a modification of Gumbel's first asymptote gives a 
perfect description of the data over about 80 orders of 
magnitude and the agreement appears to continue below 
the critical temperature. Although the detailed relation- 
ship between extreme order statistics and the EAI model 
remains to be understood, it is certainly quite rare that a 
physical formula has been tested over such a large range. 

The statistics of extremes was pioneered by Frechet, 
Fisher and Tippert, von Mises, and Gumbel. A stan- 
dard result 1^,^, due to Fisher and Tippert, Kawata, and 
Smirnow, is the universal distribution of the first, second, 
third,. . . smallest of a set of N independent identically 
distributed random numbers. For an appropriate, expo- 
nential decay of the random number distribution their 
probability densities are given by 



fa{x) = Ca exp [a (a 



)] 



(1) 



in the limit of large N. Here x is a scaling variable, which 
shifts the maximum value of the probability density to 
zero, and C'a — a°/r(a) normalizes the integral over 
faix) to one. In Gumbel's book Eq. (|l|) is called the 
first asymptote, as it holds for the asymptotic extreme 
order statistics of the first of altogether three different 
universality classes of random number distributions. The 
exponent a takes the values a — 1, 2, 3, . . ., correspond- 
ing, respectively, to the first, second, third,. .. smallest 
random number of the set. This holds independently of 
the details of the original random number distribution, 
as long as one stays within the first universality class. 
In the last years a non-integer value of the exponent a 
received also some attention. For the probability den- 
sity of the magnetization of the 2D XY model Bramwell 
et al. §,|lll derived a = 7r/2 in the spin wave approxi- 
mation and conjectured that this exponent describes, at 
least approximately, probability densities of a large class 
of correlated systems. 

For disordered systems Bouchaud and Mezard 0| 
noted that a relationship to extreme order statistics is 
intuitively quite obvious. Namely, at low temperatures 
a disordered system will preferentially occupy its low 
energy states, which are random variables due to the 
quenched exchange interactions of the system. Their 
investigation shows that Gumbel's first asymptote with 
a = I corresponds to one step replica symmetry break- 
ing, and their conjecture of a relationship between ex- 
treme order statistics and disordered systems is certainly 
far more general. This, and the possible description of a 
broad range of critical phenomena by the a = tt/2 modi- 
fication of Gumbel's first asymptote, has motivated us to 
analyze the overlap probability density of the EAI model 
at and below the critical point with respect to the large 
X fall-off behavior of Eq. (|l|) . 

The energy function of the J — ±1 EAI spin-glass 
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model is given by ||^ 



E ^ ~ ^ Jik SiSk , 
iik) 



(2) 



where the = ±1 are the spins of the system and the 
sum is over the nearest-neighbor pairs of a cubic lat- 
tice with periodic boundary conditions. The coupling 
constants Jik are quenched random variables, which take 
on the values ±1, with equal probabilities. A set of cou- 
pling constants defines a realization J = [Jik] of the 
system. The two replica overlap (Parisi order parame- 
ter) is defined by 



where the s"^^ and sY' are the spins of two copies 
(replica) of the realization J and the sum is over all sites. 
The overlap probability density is given by the average 
over the probability densities P'-^ {q) of all realizations 



1=1 



,(2) 
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where iVj- is the number of realizations used and L is the 
lattice size. There is a long history of MC studies of this 
model [^-19 1, which have led to a wealth of informa- 
tion. Here we introduce only two results: 1. The model 
has a freezing transition at a finite temperature, which 
according to the most recent estimates [^,^ is consis- 
tent with Tc = 1.14, where the temperature is given in 
dimcnsionless units such that the Boltzmann constant is 
T4|-p6 18 1 reported results which were 
~ (KT) type fine 



one. 2. References 

consistent with a Kosterlitz-Thouless 
of critical points below Tc, quite similar to the 2D XY 
model. In our context this is of interest in view of the de- 
scription of this model by Eq. (|^) with a = 7r/2 Note, 
however, that one of the most recent EAI investigations 
pof claims to rule out the KT scenario. 

At T = 1.14 we generated 8192 realizations for L = 
4, 6, 8, 512 realizations for L = 12 and 256 realization 
for L = 16. Figure |l| shows our normalized P((?) prob- 
ability densities. Due to the MC method [10| used, the 
error bars of neighboring entries are strongly correlated. 
This results in smooth curves of varying thickness, which 
represents the error. The pearl of these data are the 
tails of the distributions, which (for L = 16) are accu- 
rate down to 10^^®". This is achieved by simulating the 
system in a statistical ensemble for which the distribu- 
tion of g- values is approximately flat, instead of using the 
Gibbs canonical ensemble. After the simulation, results 
for the Gibbs ensemble are obtained through an exact re- 
weighting procedure. In this way computer simulations 
allow to probe much easier into the extremes of materi- 
als than real experiments. Alongside with our data at 
the critical point, we analyze our data from our simula- 
tions [[l7|Jl8| at T = 1, below the critical point. In that 



Q. 




FIG. 1. Overlap probability densities Pl{<1) versus q for the 
EAI model on lattices at the critical temperature. 



case we generated 8192 realizations for i = 4, 6, 8 and 
640 realizations for L = 12. In the tails the data (for 
L — 12) at T = 1 are accurate down to 10"^'^. 

We first ask the question whether, up to finite-size cor- 
rections, the probability densities depicted in Fig. || scale. 
A method to investigate this is to plot aL PL{q) versus 
(m ^ where is the mean value of q with respect 

to the distribution PL{q) and aL is its standard devia- 
tion (here (Jl = because the Phiq) are even functions). 
Visual inspection shows that the data scale indeed and 
we proceed to fit the standard deviations to the two pa- 
rameter form (Tl = ci L^^^'^ to obtain 



(3 



= 0.312(4), g = 0.32 for r= 1.14, and (5) 



= 0.230(4), Q = 0.99 for T = 1 



(6) 



Here the numbers in parenthesis denote error bars with 
respect to the last digits and Q is the goodness of fit. For 
T = 1.14 we plot in Fig. | P'{q') = PL{q)/L'^''' versus 
q' = L^/'^ q and see that the four probability densities 
collapse onto a single curve. To enlarge the scale, we re- 
strict ourselves to the g > range. The error bars of the 
lines in Fig.B are up to multiplicative factors the error 
bars of Fig. nT Not to obscure the agreement, we include 
only one representative error bar for each lattices size, 
L = 16, 12, . . . from right to left in the Fig. | (for L < 8 
they are barely visible on the scale of the figure) . For our 
data at T = 1 a similar analysis is already given in . 
The small discrepancy in the estimates of the critical ex- 
ponent P/iy (0.255 in instead of 0.230) is due to using 
different methods of data analysis. Note that the error 
bars of the P/v estimates (||) and (^) refiect only the 
fiuctuations of our two parameter fit and additional (sys- 
tematic) errors are expected from corrections to scaling. 

Our aim is to relate the probability distribution of 
Fig. H to the first asymptote of extreme order statistics. 
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FIG. 2. Rescaled overlap probability densities for the EAI 
model on lattices at the critical temperature. 



Q. 




-50 
-100 
-150 
-200 
-250 
-300 
-350 
-400 



— ^^^^ — '"^..--^ 






N;x T=1.14 


- 

L=4 


\ 

\ 


L=6 


\ 


L=8 




L=12 




L=16 




Eq.(11) 









1.5 1.6 1.7 1.8 1.9 2 

q' 



2.1 2.2 2.3 2.4 



FIG. 3. Tails of the rescaled overlap probability densities 
of Fig. I ln[P'(g')] versus q' . The error bars are smaller or 
of the order of the thickness of the lines. 



Eq. In the neighborhood of x = 



x-e' = -l--x^ + 0{x^) 



(7) 



holds. To get the position of the maximum of P'{q') 
right, we have to choose 



x = h{q' - q'^^J , 



(8) 



where g^^ax the q' argument for which the probability 
density P'{q') takes on its maximum value. We then 
have to find an exponent a to reproduce the data for 
x^b{q' - q'^^^ > 0. For a; = 6 {q' - q'^^^^) < 0, however, 



the behavior cannot be quite correct. The reason is 
that the a; < asymptotic behavior 



exp [ a a; ] = exp [ a 6 (g' - g,'„ax) ] 



(9) 



predicts on a logarithmic scale a constant slope a with 
decreasing x, while for the data of Fig. || the slope levels 
off and at a; = —bq'^^^^ we have 



d \n[P'{q') 



dq' 



= 0, 



(10) 



q'=0 



what is impossible with (^. A simple solution is to re- 
place the first X on the right-hand side of Eq. (|l]) by 
c tanh(a;/c), where c > is a constant. For small x 
the Taylor expansion still holds, while for large |x| 
the hyperbolic tangent function c tanh(x/c) approaches 
quickly ±c (note that in the limit c —>■ oo the original 
form (|^) is recovered). For a; ^ — oo (practically already 
at q' — 0) the thus modified Gumbel distribution be- 
comes constant. Therefore, the symmetric expression for 
P'{q') is obtained by multiplying the above construction 
with its reflection about the q' = axis. 



P'{q') = C X 
exp < a 



(11) 



ctanh(+J(g'-g;,,^^l-.+''(^'-«". 



exp < a 



c tanh ( —-((?' + q[a 



Of course, the important large x behavior of Eq. (^ is 
not at all affected by our manipulations. 

The calculation of the parameters a, b, c and C is done 
by using the logarithm of Eq. ([l^). The constants are 
essentially determined by the following parts of the dis- 
tribution: C by the height of the peak at q' — q'^^^^^ Wmax 
is now off the maximum location by a tiny shift which 
can be neglected; in the fits we used g^ax = 1.135972 for 
T = 1.14 and q[^^^^ = 1.115056 for T = 1), a and 6 by 
the q' > (/inax tails of the distribution and c by the value 
at q' = 0. This allows to iterate from a reasonably close 
starting values to our final estimates 



a = 0.448 (40) for T = 1.14, 
a = 0.446 (37) for T = 1 . 



and 



(12) 
(13) 



The non-universal coefficients are b = 5.35 (11), c — 
3.37 (41), C = 7.39 (86) for T = 1.14, and b = 8.23 (17), 
c = 4.48 (43), C = 16.8 (2.2) for T = 1. The error bars 
rely on a jackknife analysis. 

For T — 1.14 our best fit to Eq. ( |ll| ) is already included 
in Fig. ||. In Fig. || we follow the tails of our distributions 
by plotting ln[P'((7')] versus q' for q' > 1.5. Besides T = 
1.14, the results for T = 1 are also included in this figure. 
On the scale of Fig. || the error bars are smaller than the 
thickness of the lines, namely in the range 0.01 to 3.1 
(exceeding 1.0 only for the largest q' values of the L = 12 
and 16 lattices). The reason is that our sampling is flat in 
q and the huge range of ln[P(g)] comes from subtracting 
logarithms of (numerically exact) weight factors. 
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Figure ^ exhibits the finite-size effect that for q close 
to one the smaller lattices undershoot the larger ones. 
It is quite clear that something like this has to hap- 
pen, because the data from each lattice size terminates 
at g = 1, whereas Eq. (^) has no corresponding singu- 
larity. When calculating our fit parameters, we take this 
into account by restraining our use of data to q' < 2, 
{\n[P{q')] > -43.4), for the T = 1.14, L = 16 lattice 
and to q' < 1.62, (ln[P(g')] > -25.6), for the T = 1, 
L = 12 lattice. The agreement of our fits with those 
data stretches then over considerably larger ranges. Dis- 
crepancies of the L = 16 data with the fit begin only 
around ln[P(g')] = —200, and even this is only visible on 
a major enlargement of the scale of Fig. ||. Discrepancies 
of the r = 1, L = 12 lattice with the fit are encountered 
around ln[P(g')] = -35. However, the T = 1.14, L = 12 
data deviate already around ln[P((7')] = — 10 from the 
L = 16 data. This appears possible, because corrections 
to the scaling factor are not traced by the accuracy 

of our data (in particular L = 16 has low statistics due 
to computer time limitations). Therefore, it is not en- 
tirely clear whether the large range which we find for the 
agreement of the fit with our L = 16 lattice is to some 
extent a statistical accident. Taking it at face value, we 
have the remarkable range of 200/ln(10) ~ 87 orders of 
magnitude. An actually smaller range would still be large 
enough to give us confidence that we are dealing with a 
true effect. 

Our coefficient a is far off from the 2D XY coefficient 
of Bramwell et al. a = 7r/2. This means that the EAI 
and the 2D XY models are certainly in quite different 
universality classes of extreme order statistics. However, 
the fact that both distributions can be described by it at 
all might help to explain the observed similarities. Our 
temperature P = 1 is below the critical Tc, but with our 
lattice sizes it appears impossible to resolve the question, 
whether the here reported behavior reflects the existence 
of a critical line below Tc or just the closeness of T = 1 to 
Tc- Before comparing to extreme order statistics, we fl^ j 
tried to fit the q > Qmax tails of our distributions to the 
theoretical predictions which have been made by Parisi 
and collaborators |^-^ based on the replica mean-field 
approach. None of these fits was particularly good and 
even when pushing the adjustment of free parameters to 
their limits only small parts of the tails of our distribu- 
tions could be covered. 

In summary, we have presented strong numerical ev- 
idence that the Parisi overlap distribution of the EAI 
model can be described by Eq. (|l^). The detailed re- 
lationship between the EAI model and extreme order 
statistics remains to be investigated and it is certainly 
a challenge to extend the work of Bouchaud and Mezard 
to the more involved scenarios of the replica theory. 
On the other hand, it could be that replica symmetry 
breaking is not the driving mechanism of the EAI model 
phase transition and that our observations are rooted in 
general relations Q between correlated systems and ex- 
treme order statistics. 
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